import ij.IJ;
import ij.ImagePlus;
import ij.ImageStack;
import ij.gui.GenericDialog;
import ij.plugin.filter.PlugInFilter;
import ij.process.ByteProcessor;
import ij.process.FloatProcessor;
import ij.process.ImageProcessor;


public class A07_Gaussian_Average implements PlugInFilter{

	ImagePlus img1;
	double TAU = 50;
	double ALPHA = 0.05;
	double sigma0 = 10.0;
	double S = 4;
	
	@Override
	public int setup(String arg0, ImagePlus img) {
		img1 = img;
		return DOES_8G + STACK_REQUIRED + NO_CHANGES;
	}
	
	@Override
	public void run(ImageProcessor ignore) {
		if(img1 == null)
		{
			IJ.log("Kein Stack");
		}
		
		if (!getUserInput()) {
			return;
		}
		
		// Get stack1 from the source image:
		ImageStack stack1 = img1.getImageStack();
		
		ImageStack µk = new ImageStack(stack1.getWidth(), stack1.getHeight());
		ImageStack sk2 = new ImageStack(stack1.getWidth(), stack1.getHeight());
		ImageStack Mk = new ImageStack(stack1.getWidth(), stack1.getHeight());
		
		//Making BackgroundModels
		for(int i = 1; i <= stack1.getSize(); i++)
		{
			if(i == 1)
			{
				ImageProcessor ip = stack1.getProcessor(i);
				µk.addSlice(ip);
			} else {
				ImageProcessor ik = stack1.getProcessor(i);
				ImageProcessor µ = µk.getProcessor((i-1));
				ImageProcessor ip = new ByteProcessor(stack1.getWidth(), stack1.getHeight());
				
				for(int u = 0; u < ip.getWidth(); u++)
				{
					for(int v = 0; v < ip.getHeight(); v++)
					{
						double value = ALPHA * ik.getPixel(u, v) + (1- ALPHA) * µ.getPixel(u, v);
						ip.putPixel(u, v, (int)value);
					}
				}
				µk.addSlice(ip);
			}
		}
		
		(new ImagePlus("Gaussian average µk", µk)).show();
		
		//Sigma k quadrat
		for(int i = 1; i < stack1.getSize(); i++)
		{
			ImageProcessor ik = stack1.getProcessor(i);
			ImageProcessor µ = µk.getProcessor(i);
			if(i == 1)
			{
				ImageProcessor ip = new ByteProcessor(stack1.getWidth(), stack1.getHeight());
				for(int u = 0; u < ip.getWidth(); u++)
				{
					for(int v = 0; v < ip.getHeight(); v++)
					{
						ip.putPixelValue(u, v, sigma0 * sigma0);
					}
				}
				sk2.addSlice(ip);
			} else
			{
				ImageProcessor ip = new ByteProcessor(stack1.getWidth(), stack1.getHeight());
				ImageProcessor sigmaM1 = sk2.getProcessor(i-1);
				for(int u = 0; u < ip.getWidth(); u++)
				{
					for(int v = 0; v < ip.getHeight(); v++)
					{
						double value = ALPHA * Math.pow((ik.getPixel(u, v) - µ.getPixel(u, v)), 2) + (1 - ALPHA) * sigmaM1.getPixel(u, v);
						ip.putPixelValue(u, v, value);
					}
				}
				sk2.addSlice(ip);
			}
		}
		
		(new ImagePlus("Gaussian average sk2", sk2)).show();
		//Mk		
		for(int i = 1; i < µk.getSize(); i++)
		{
			ImageProcessor ip = new ByteProcessor(stack1.getWidth(), stack1.getHeight());
			ImageProcessor ik = stack1.getProcessor(i);
			ImageProcessor µ = µk.getProcessor(i);
			ImageProcessor sk = sk2.getProcessor(i);
			
			for(int u = 0; u < ik.getWidth(); u++)
			{
				for(int v = 0; v < ik.getHeight(); v++)
				{
					if(Math.abs(ik.getPixel(u, v) - µ.getPixel(u, v)) > S * Math.sqrt(sk.getPixelValue(u, v)))
					{
						ip.putPixel(u, v, 255);
					} else {
						ip.putPixel(u, v, 0);
					}
				}
			}
			Mk.addSlice(ip);
		}
		
		// Display stack2:
		(new ImagePlus("Gaussian average", Mk)).show();
	}

	private boolean getUserInput() {
		GenericDialog gd = new GenericDialog("Select threshold");
		gd.addSlider("Tau", 0, 255, TAU);
		gd.addSlider("ALPHA", 0.01, 0.99, ALPHA);
		gd.addSlider("Sigma0", 0.01, 100, sigma0);
		gd.addSlider("S", 0.01, 10.00, S);
		gd.showDialog();
		if(gd.wasCanceled())
		{
			return false;
		}
		TAU = gd.getNextNumber();
		ALPHA = gd.getNextNumber();
		sigma0 = gd.getNextNumber();
		S = gd.getNextNumber();
		return true;
	}

	

}
